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Abstract. We present a numerical code for multi-component simulation of the galactic 
evolution. Our code includes the following parts: iV-body is used to evolve dark matter, stellar 
dynamics and dust grains, gas dynamics is based on TVD-MUSCL scheme with the extra 
modules for thermal processes, star formation, magnetic fields, chemical kinetics and multi¬ 
species advection. We describe our code in brief, but we give more details for the magneto-gas 
dynamics. We present several tests for our code and show that our code have passed the tests 
with a reasonable accuracy. Our code is parallelized using the MPI library. We apply our code 
to study the large scale dynamics of galactic discs. 


1. Introduction 

An explosive growth of observational data on dynamics and chemical composition of various 
components (subsystems) of galaxies is seen in last decade. To interpret the huge data cubes we 
need reasonable models included more physical processes like magnetic fields, star formation, 
feedback, chemical kinetics, dust grain dynamics and so on. Spatial scales of the processes in a 
galaxy vary from sub-parsecs for turbulent structures in the warm interstellar medium to several 
kiloparsec size for galactic disc and winds. Formation of molecules and distribution of dust grains 
strongly affect on the efficiency of star formation, which, in turn, changes the thermal state of a 
gas through stellar winds and supernova explosions [1]. Then, to study the interaction between 
various components numerically someone needs a code with a lot of physical modules. 

In recent years several open source codes were developed for astrophysical applications, such 
as ZEU^ RAMSE^ GADGET]^ Athens]^ FLASH|^and many others. These codes are general- 
purpose ones, so that each of them has own advantages and disadvantages in the application 

^ http://www.astro.princeton.edu/ jstone/zeus.html 

^ http://irfu.cea.fr/Phocea/Vie_des_labos/Ast/ast_sstechnique.php?id_ast=904 
^ http://www.mpa-garching.mpg.de/gadget/ 

^ https://trac.princeton.edu/Athena/ 

^ http://flash.uchicago.edu/site/ 


to the galactic evolution [21E1I11IS]. However galaxies are multi-component systems, and its 
simulation requires a robust model included special physical processes. So the main goal of this 
project is to develop the numerical code for simulation of the evolution of disc galaxies including 
generation of spiral structure, physics of interstellar medium, formation of clouds and stars ©C]. 

This paper is organized as follows. In Section 2 the main numerical methods are described. 
Section 3 presents the results of tests. Section 4 describes the application of our code to the 
evolution of stellar-gaseous galactic disc. In Section 5 we summarize our results. 


2. Description of the code 
2.1. Gas dynamics 

The gas dynamics equations for an ideal magnetized gas can be written in the conserved 
variables U in Cartesian coordinates 


U — [p, pVx^) P'^y-) P'^Z') 

Then, the conservation laws can be written in a compact form 

dU dF dG dU 

dt dx dy dz 
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where F, G, and H are vectors of fluxes in x, p, and 2 ;-directions, respectively and S is the 
source vector. The components of these vectors are 
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For pure gas dynamics we 


In multi-dimensional magneto-gas dynamics to obtain zero-value for the magnetic field 
divergence we apply the Constrained Transport technique for magnetic field transport through 
computational domain BM- In this approach magnetic field strength is defined at faces of a 
cell, while other gas dynamical variables are defined at the center of a cell. 

We solve the equations (§ in Cartesian coordinates [x^y^z). The center of each cell is 
located at the position (x^, pj, Zk)- For each cell the faces being normal to the x direction have 
coordinates x^ii/2 ^md sizes 5x, 5y and 5z. The gas dynamic variables (density, momentum, 
energy) are volume-averaged and defined at the center of a cell. For example, the density value 
can be written as follow 
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While the magnetic field components are surface-averaged and found at faces of a cell: 
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To calculate the momentum and energy fluxes we apply the second-order cell-centered averages 
for magnetic field, e.g. for Bx\ 

'^{^x,i-\-l/2J,k T l/2,j,fc) • (6) 

We use a finite-volume discretization for gas dynamics equations and finite-area discretization 
for magnetic field transport equations. In this case gas-dynamical variables at time t + 6t can 
be found as follow 

TTn+l _ TT^ _ Tp^+1/2 _ -p,n-|-l/2 _ ^n-\-l/2 _ ^n-\-l/2 _-pj-n-|-l/2 _ -pj-n-|-l/2 

^ij,k — ^iJ:k i^l/2J,k ^i-l/2J,k\ gy [ z,j+l/2,fc ^ij-l/2,k\ [ z,j,fc+l/2 ^iJ,k-l/2\ ’ 

(7) 

where the vectors ^rj+ 1/2 fc’ '^^]^k+i/2 time- and area-averaged fluxes through 

the X, y and z-faces of a cell, correspondingly. 

Using the Stokes’ law for the magnetic field transport equation we obtain time-evolution of 
magnetic field equation: 

oTi+l _ TDn _ pn-\-l/2 _ pn-\-l/2 , pn-\-l/2 _ pn-\-l/2 
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( 8 ) 

where £x, £y and £z are the components of electric field (or electro-magnetic force): 

^ = -V X B . (9) 

The approximation of the electro-magnetic force can be written as follow: 

^x,j-l/2,A:-l/2 = ^ ^x,j-l/2,fc + ^x,j-l/2,k+l + ^x,j,/c-l/2 + £x,j+l,k-l/2 + 
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where the derivative of £x for each grid cell face is computed by selecting the upwind direction 
according to the contact mode, namely, 
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Note that in 3D the expressions similar to the above are required to convert the x- and y- 
components of the electric field to the appropriate cell corners. These expressions might be 
obtained directly cyclic permutation of the (x^y^z) and 

Godunov-type methods are widely used for solving of hyperbolic equations. An idea is to 
use the Riemann solution for the decay of an arbitrary discontinuity. In this case it is assumed 
that the solution might be non-smooth and discontinuities might be in every computational cell. 
That is why, this kind of numerical schemes is a universal tool for simulation of gasdynamical 
problems of evolution shocks and contact discontinuities. The main problem of Godynov-type 
methods is a necessity of solving system of nonlinear equations. Usually this leads to iterative 
procedure, which requires additional computational resources. To resolve this issue there are 
a lot of approximate solvers of Riemann problem. Our code includes two types of solvers for 
pure gas dynamics: Harten-Lax-van Leer (HLL) or Harten-Lax-van Leer-Contact (HLLC), and 



one for magneto-gas dynamics, namely, Harten-Lax-van Leer-Discontinuities (HLLD). These 
methods are described by m in detail. 

An exact or approximate Riemann solver requires values for the conserved variables at 
left and right interfaces of each cell: Ul — Ui — 1/2, Ur — Ui — 1/2. These values can be 
reconstructed from cell centered values with second-order or third-order accuracy, on a choice. 
The interpolation of the values is carried out for primitive variables w = {p, P, v,B} as follow 


dwi = Wi^i - Wi , ( 12 ) 

dWi = K2io[immod{dwi-i^ dwi) -h Piminmod((ii(;^, bdwi-i ), (13) 

dWi^i = K2io[immod{dwi^i^ dwi) -h Piminmod((ii(;^, bdwi^i ), (14) 

and then wr = Wi + dWi wr = Wi^i — dWi^i , where Ki = ^(1 — a^), K 2 = ^(1 + a^). 


b — {3 — r) /{I — r) and minmod is a limiter function [TT]. The parameter r can be equal 
to —1, 0, or 1/3, which represent second-order fully upwind, second-order upwind biased, and 
third-order upwind biased cases, respectively. Below we use a^ = 1/3. 

2.2. N-body dynamics 

We use the ’’particle-in-cell” method for the dynamics of both dark matter halo (’’live” halo) 
and stellar discs. To obtain the second-order accuracy on time we apply the ’’flip-flop” integrator. 
The main advantage of this approach is that it requires only one calculation of forces per time 
step for each particle. More details about our Wbody solver for galaxy dynamics can be found 
in [12]. 

2.3. Poisson solver 

To take into account self-gravity of a gas we solve the Poisson equation with the gas dynamics 
conservation laws (§: 

(x, y, z) = Ai^Gp{x, y,z). (15) 

In our code the Poisson equation can be solved using three different methods: FFT-based, 
TreeCode [I3| and over-relaxation Gauss-Seidel method. To calculate self-gravity in N- 
body/gasdynamical problems we use regular Cartesian grid. For interpolation of particle density 
to a mesh we use the second order finite-volume interpolation, which conserves the total mass 
of the system with a good accuracy. 

3. Test problems 

Below we present several tests for our gasdynamical code in ID, 2D and 3D. 

3.1. The Sod shock tube test 

This is the most simple test demonstrated the formation of shock wave, contact discontinuity 
and rarefaction wave. The initial distribution of gasdynamical values is taken as usual: 
p, Vx^ P = {1, 0,1} at X < 0; p, Vx^ P = {0.125, 0,0.1} at x > 0. Initially the velocity is equal to 
zero. Figure [^presents the result of numerical simulation for different resolution and the exact 
solution. One can see a good coincidence of the numerical solution with the exact one even for 
low resolution. 






Figure 1. The Sod shock tube. The exact solution is shown 
by black line, the other lines represent the distribution at 
t = 0.5 with grid resolution N — 200 (green), N = 500 (red), 
N = 1500 (blue). One can find five regions in the solution: 
the regions 1 and 5 correspond to the initial unperturbed 
state of a gas, the rarefaction wave is found in the region 2, 
the regions 3 and 4 are separated by contact discontinuity, 
and the shock wave is located between regions 4 and 5. 


Figure 2. The Shu-Osher 
test. The initial density dis¬ 
tribution is shown by black 
line, the other lines represent 
the distribution at t = 0.18 
with grid resolution N = 200 
(green), N — 500 (red), N = 
1500 (blue). 


3.2. The Shu-Osher problem 

The Shu-Osher problem tests a shock-capturing scheme’s ability to resolve small-scale 
flow features. It shows the numerical viscosity of the method. The initial distribution is: 
p = 3.857143, p = 10.33333, u = 2.629369 at x < 0.125 and p = 1 -\- 0.2 sin( 167 rx), p = 1, = 0 

otherwise, 7 equals 1.4. 

Figure shows the initial density distribution (black line). We study the dependence on 
spatial resolution: the case for N = 200 is exactly smooth, whereas oscillations appear for 
N = 500, for further increase of resolution (N = 1500) our result becomes closer to the exact 
solution M- 

3.3. The bow shook simulation 

A bow shock usually forms due to supersonic motion of star (or planet) through interstellar 
medium |15j or a galaxy through the intercluster medium 1E]> SO that this is one of the 
most important astrophysical problem. To test our code we simulate a supersonic gas flow 
through potential well. The initial distribution of the parameters are homogeneous: p = 1.0, 
P = 1.0, = 1.0, and 7 equals 1.4. The external potential well is set in the form 

4^(x) = — 4 ^ 0 exp(—(x/xo)^), where 4^o = 5 and xq — 0.1. The boundary conditions are free 
at the right and the steady inflow at the left. 

At t = 0.051 the shock wave is formed on the far edge of the well (towards to the flow) due 
to supersonic falling of the gas to the potential well (red line). Note that this configuration is 
unstable. So the shock wave moves upstream through the potential well and after crossing the 
minimum of the potential it becomes steady-state at front edge of the potential well (blue line). 

5 . 4 . The Brio-Wu problem 

The Brio-Wu test is one dimensional magneto-gas dynamics problem. The solution of this 
test consists of the fast rarefaction wave, that moves to the left, the intermediate shock wave, the 
slow rarefaction wave, the contact discontinuity, slow shock wave and another fast rarefaction 
wave, that moves to the right (see figure Q. The intermediate shock and slow rarefaction waves 

























Figure 3. The bow shock test. The density distribution is at t = 0.051 (red), 0.226 (green) 
and 0.375 (blue). The potential well T(x)/To is shown by black solid line. 




Figure 4. The Brio-Wu problem. Our result is depicted by red dots, the result obtained in the 
Athena code is shown by blue line. 


form a structure called the compound wave m- The initial distribution for this test is: p = 1, 
P=l, Sy = latx<0, p = 0.125, P = 0.1, By — —1 for x >0. The longitudinal component 
of magnetic field = 0.75 is constant over the grid. 

This test (and the next one, the Ryu-Jones problem) is compared with the results obtained by 
the Athena code Hi. For both simulations we set the same spatial resolution N = 200 (figure Q. 
One can see a good agreement between the results, but the artificial oscillations can be found 
for the longitudinal velocity component for both solutions. 

3.5. The Ryu-Jones problem 

The Ryu-Jones problem tests the rotation of the magnetic field components. The solution 
consists of two fast shock waves with velocities 1.22 and 1.28 Mach numbers, two slow shock 








































Figure 5. The Ryu-Jones problem. Our result is depicted by red dots, the result obtained in 
the Athena code is shown by blue line. 



Figure 6. The rotor problem. The maps of density, pressure, magnetic pressure and Mach 
number are shown from left to right correspondingly. 


waves with 1.09 and 1.07 Mach numbers, two rotational and one contact discontinuities m- 
The initial distribution is: p — 1.08, P — 0.95, Vx = 1.2, Vy = 0.01, Vz = 0.5, Bx = 4/\/4^, 
By = 3.6/V^, Bz = at X < 0, and p = P = Vx — Vy — Vz = Bx — 4/\/4^, 

By — 4/V^, Bz — at X > 0. One can see a good agreement between results obtained 

both in our code and in the Athena code. However the peaks of Vy and By in our simulation are 
smoother than these obtained in the Athena code m- 

3.6. The rotor problem 

This test demonstrates fast rotation of the cylinder in the non-moving medium with 
homogeneous magnetic field. Initially there is a disc in the center of computational domain: 
p = 10, Vx = -V[)y/r^ and Vy = v^x/r^ at r < ro = 0.1; p ^ 1 + 9/(r), Vx = -VQf{r)y/r^, 
= '^o/(^)^/'^o at ro < r < ri = 0.2, and p = 1 at r > ri; P = 1 and Bx — 5/V^^ are 
constant in the whole computational domain. Here the function /(r) = (ri — r)/(ri — ro) is 
the linear interpolation of the gas density and velocity. This configuration is strongly unstable, 
because the centrifugal forces are not balanced. A rotating gas is redistributed gradually within 
the computational domain capturing the stationary gas. Under these conditions the magnetic 
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Figure 7. The Sedov-Taylor blast-wave. The 2D density map for pure gas dynamics (left 
panel) and the radial pressure profile at t = 0.1 for our simulations (red dots in left middle 
panel) compared with the exact solution (blue line). The 2D density and magnetic pressure 
maps at t = 0.1 for magneto-gas dynamics (two right panels). 


field keeps the rotating material in the flattened form (Figure]^. 

3.7. The Sedov-Taylor blast-wave 

One of the most well-known self-similar solution describes a strong shock wave originated 
from the point explosion. For zeroed magnetic field a numerical solution can be compared 
with the exact solution. In 2D we simulate a strong explosion with the following parameters: 
p = 1, P = 10“^ are set in the whole computational domain, P = 7r^rQ/(7 — 1) at r < tq. We 
set r — 0.01 and take a computational domain [—0.5,0.5] x [—0.75,0.75] and number of cells 
512 X 768. Figure]^ shows the 2D density map (left panel) and the radial pressure profile for our 
simulations (red dots in left middle panel) compared with the exact solution (blue line). One 
can see a good agreement between numerical data and analytic curve. 

We consider a blast wave in a magnetized medium. At t = 0 we set Bx — By — l/\/8^. 
Despite of the symmetric initial distribution the shock front expands larger along the lines of 
the magnetic field. 

3.8. The Kelvin-Helmholtz instability 

Shear instabilities are very common in astrophysical objects, so that numerical algorithms 
should reproduce such inabilities well. Here we simulate the evolution of two gas layers moving 
with different velocities and separated initially by a contact discontinuity: p — Vx — —0.5 
at \y\ < 0.25 and p = 2^ Vx — 0.5 otherwise, and P = 2.5 everywhere. We add a random 
perturbation of the velocity with the amplitude 0.001. The grid size is 512 x 512. 

Figure shows the evolution of gas density for pure gas dynamics (upper row of panels) 
and magneto-gas dynamics (bottom row of panels). One can note the turbulent flow formation 
due to shear instability at early times, which produces large scale vortexes at further nonlinear 
stage. The presence of magnetic field strongly changes the picture: a small scale perturbations 
disappear for initial longitudinal magnetic field Bx = 0.5 (see bottom line in Figure]^, whereas 
a large scale oscillation of a gas along magnetic field can be seen. 

3.9. Spherieal eollapse of gravitating gas 

To test a gravity solver we use a standard 3D cosmological problem - the collapse of a 
gaseous sphere with initial density profile p{r) — 1/r. We set the cubic computational domain 
{x^y^z) G [—1,1] with cell number 100^. Figure ^ shows the evolution of the density profile. 















Figure 8. The Kelvin-Helmholtz instability. The gas density distributions at t = 0,0.1,4,12 
(from left to right) in the pure gas dynamics (top row of panels) and the same, but in the 
presence of magnetic field (bottom row of panels). 


One can see that dense core forms due to homogeneous collapse and strong shock wave moves 
outwards (see right panel of the Figure). 


3.10. Thermodynamical module 

To study the evolution of the interstellar medium in galaxies we develop a module of thermal 
processes. In this module we can switch between tabulated cooling/heating rates and calculation 
of rates for main chemical species in the interstellar medium. For the former we can use tables 
of cooling rates in the temperature range 10 K< T < 10^ K calculated by [2Ql [2T1 [22]. In the 
latter we calculate the cooling and heating rates produced by specific emission processes and in 
the temperature range T < 2 x 10^ K, that is of a great importance for studying thermal state of 
the interstellar medium (see |23| for instance). We can investigate the thermal evolution of the 
interstellar medium for different abundances of heavy elements (metallicities). Our thermal block 
includes following cooling processes: cooling due to recombination and collisional excitation and 
free-free emission of hydrogen m, molecular hydrogen cooling |25| , cooling in the fine structure 
and metastable transitions of carbon, oxygen and silicon |26| . energy transfer in collisions with 
the dust particles m and recombination cooling on the dust[28|. The heating rate takes into 
account photoelectric heating on the dust particles 123 EH, heating due to H 2 formation on dust 
and the H 2 photodissociation |29] and the ionization heating by cosmic rays m- For multi-level 
atoms the cooling rates are obtained from the level population equation assuming the optically 
thin steady-state regime (see e.g. 


Figure 10 presents cooling rates in the low temperature (T < 2 x 10^ K) range adopted for 


our model of the interstellar medium. 















































Figure 9. The density profile for t = 0 and 
t = 1 is depicted by cyan and black lines, 
respectively, the pressure profile for t = 1 is 
shown by red line (left panel). The velocity 
profile for t = 1 is shown at right panel. 
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Figure 10. Cooling functions, A/n|^, for solar 
metallicity. The total cooling rate is depicted 
by black line. 


3.11. H 2 kinetics 

Molecular hydrogen (H 2 ) regulates star formation in galaxies [23], then to study the evolution 
of the interstellar medium we need to take into consideration the H 2 kinetics. Because of 
significant difference between dynamic and chemical timescales the H 2 kinetics is nonequilibrium 
and to get a source term in equation (§ we solve a system of ordinary differential equations 
for the following chemical species: H, H^, H 2 . Molecular hydrogen is formed on the surface 
of dust grains and dissociated by ultraviolet Lyman-Werner photons and cosmic rays. Because 
of strong absorption and scattering of ultraviolet photons in the neutral hydrogen and dust to 
mimic radiative transfer and H 2 self shielding we use a simplified approach introduced by m 
for calculation of neutral and molecular hydrogen column densities. 

3.12. Star formation implementation 

Stellar feedback effects can significantly influence on the galactic evolution. However, there 
is no consensus about the implementation of the star formation into gas dynamics, that is 
clearly demonstrated in [32]. Their results of the galaxy evolution simulations within the 
ACDM paradigm strongly depend on the star formation and feedback recipes: stellar mass, 
size, morphology and gas content of the galaxy at z = 0 vary significantly due to the different 
implementations of star formation and feedback. Despite this problem we include the star 
formation effects using an intuitive approach, that partly based on the well-known method 
offered in [33] . 

To form a star or in general stellar particle we need to find cells in a computational domain, 
which satisfy to our conditions for star formation. These conditions can be more or less 
sophisticated and sometimes may have intricate nature. Usually the following criteria for the 
stellar particle creation in a given cell are considered: the surface density should be greater 
than the threshold for star formation and simultaneously the temperature should be less 
than some critical level Tf. In the simulations presented below we have assumed > 10^ cm“^ 
and Tf < 10^ K. If these conditions are realised in a given cell, then we create a test stellar 
particle. The Salpeter initial mass function is assumed for each particle, which usually consists 
of 10^ — 10^ stars. The initial velocity of this test particle is taken to be the same as its host 
cell. The kinematics of test particles are followed by the second-order integrator mentioned in 
Section 2.2. We take into account several feedback effects: stellar winds from massive stars, 
radiative pressure, supernova explosions, and stellar mass loss by low-mass stars [3l]. Metals 
ejected by supernovae and lost by low-mass stellar population can strongly change chemical 














































Figure 11. The Galactic disc evolution. First three columns: stellar (top) and gaseous (bottom) 
distributions at t = 0, 300, 600 Myr. Right pictures: radial variations of the module of magnetic 
field (top) and vectors of magnetic field in the disc plane. 


composition and thermal evolution of a gas. Because of local character of enrichment process 
the re-distribution (mixing) of metals is of great importance for further star formation in 
galaxies [35l |36l [371 ISH] • 

4. Simulation of the multicomponent galactic disc 

Here we present an example of the gravitationally unstable stellar-gaseous galactic disc 
evolution. In this simulation we take into account magnetic field and radiation processes. The 
stellar disc is simulated using A^-body method with number of particles N = 10^. The number of 
cells for gas dynamics equals to 1024 x 1024 x 128. We start our simulation from the equilibrium 
state of the galactic gaseous disc within the fixed potential of the dark matter halo for Toomre’s 
parameter Qt = 1.2. A construction of initial equilibrium configuration for stellar-gaseous disc 
is described by [39] in detail. We assume that the galactic disc consists of stellar component 
and gas with temperature 10^ K in the middle plane. The magnetic field is set as the toroidal 
configuration with B — Q jiG. 

Figure [TT| presents the result of this simulation. The distribution of stellar particles (top row 
of three left panels) and gas surface density (bottom row of three left panels) are shown at time 
moments t = 0; 300; 600 Myr. One can see how the gravitational instability drives the formation 
of spiral pattern, which is clearly seen in both stellar and gaseous components. At t — 300 Myr 
it is easily found narrow and dense galactic shocks (bottom row of panels). Because of relatively 
small initial value Qt a flocculent structure is rapidly developed in the spiral arms jlO]. The 
spatial structure of a gas becomes very complex due to numerous shear flows, thermal and 
gravitational instabilities, magnetic field pressure influence as well. 

One can see the formation of the large number of dense gaseous complexes with number 
density ^ 300 — 1000 cm“^. These complexes are mostly located in galactic spiral arms. 
This result is almost independent on numerical resolution and slightly dependent on the initial 
magnetic field. The increase of numerical resolution leads to formation of clouds on smaller 



















scales. Such clouds are expected to give birth of stars. 

Figure [T^ upper right panel) shows the structure of magnetic field at galactic disc at 
t = 600 Myr. The amplitude of magnetic field decreases with radius from \B\ — 8fiG at 
4 kpc down to \B\ = 4/iG at the outer part of the disc. This is in a good agreement as with 
observations m and with cosmological simulations m- The vector map (right bottom panel in 
Figure [TT|) represents the chaotic and regular velocity components. The regular field corresponds 
to the rotation of a gas and large-scale spiral structure in the disc. 

5. Conclusion 

In this paper we have described our three-dimensional numerical code for multi-component 
simulation of the galactic evolution. This code has been mainly developed to study the evolution 
of disc galaxies taking into account generation of spiral structure, physics of interstellar medium, 
formation of clouds and stars. Our code includes the following ingredients: A^-body dynamics, 
ideal magneto-gas dynamics, self gravity of gaseous and stellar components, cooling and heating 
processes, star formation, chemical kinetics and multi-species gaseous and particle (for dust 
grains) advection. We present several tests for our code and show that our code have passed 
the tests with a resonable accuracy. It should be noted that our code is parallelized using the 
MPI library. We apply our code to study the large scale dynamics of galactic discs in context 
of formation and evolution of galactic spiral structure [39] , molecular clouds formation 031 and 
disc-to-halo interactions M- 
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